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Abstract 

We investigate smooth and sparse optimal control problems for convective FitzHugh- 
Nagumo equation with travelling wave solutions in moving excitable media. The 
cost function includes distributed space-time and terminal observations or targets. 
The state and adjoint equations are discretized in space by symmetric interior point 
Galerkin (SIPG) method and by backward Euler method in time. Several numer¬ 
ical results are presented for the control of the travelling waves. We also show 
numerically the validity of the second order optimality conditions for the local 
solutions of the sparse optimal control problem for vanishing Tikhonov regular¬ 
ization parameter. Further, we estimate the distance between the discrete control 
and associated local optima numerically by the help of the perturbation method 
and the smallest eigenvalue of the reduced Hessian. 

Keywords: FitzHugh-Nagumo equation; traveling waves; sparse controls; 
second order optimality conditions; discontinuous Galerkin method. 


1. Introduction 

Spatially and temporally varying structures occur in form of Turing patterns, 
traveling waves, fronts, periodic pulses in many physical, chemical, and biolog¬ 
ical systems. They are described mathematically in form of coupled semi-linear 
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partial differential equations (PDEs) |l30|. The FitzHugh-Nagumo (FHN) equa¬ 
tion is one of the most known generic model in physiology, describing complex 
wave phenomena in excitable or oscillatory media. The most known type of the 
FHN equation in the literature consists of a PDF with a non-monotone nonlinear 
term, known as activator equation, and an ordinary differential equation (ODE), 
known as inhibitor equation. We call such kind of activator-inhibitor system as 
classical FHN equation. Another type of the activator-inhibitor system is diffu¬ 
sive FHN equation consisting of an activator PDF, and one or two inhibitor PDEs 
II25II . Recently, the convective FHN equation has been proposed as a model for 
wave propagation in blood coagulation and bioreactor systems ifTSlfT^I^ . The 
presence of convective field leads to complex wave phenomena, like triggering 
and autonomous waves in a moving excitable media ifT^ . 

The classical FHN equation and its PDF part the Schlogl or Nagumo equations 
were investigated theoretically and numerically for the wave-type optimal control 
solutions |I71[To1[IT1|38]|. Optimal control of semi-linear parabolic equation is an 
active research field with many applications in controlling pattern formation ll^Tll 
and feedback control of the monodomain equations in cardiac electrophysiology 
ll6ll to give a few examples. 

In this paper, we investigate the numerical treatment of optimal control prob¬ 
lems governed by the convective FHN equation. The uncontrolled solutions of 
the convective FHN equation behave like travelling waves ifT^ . To control such 
a travelling waves, we use sparse control, which is a non-smooth -control cost 
in addition to the L^-control cost. When the control signals are localized in some 
regions of the space-time cylinder, sparse control provides solutions without any a 
priori knowledge of these sub-areas. Sparse optimal control was first investigated 
in II40II for linear elliptic equations and later studied in ||9l HH HU for semi-linear 
elliptic and parabolic equations. 

Here, we use symmetric interior penalty Galerkin (SIPG) method for space 
discretization. The discontinuous Galerkin (dG) methods are more stable for con¬ 
vection dominated problems than the continuous finite element methods and they 
do not require the stabilization terms like the streamline upwind/Petrov-Galerkin 
method (SUPG). The dG methods have several advantages compared to other nu¬ 
merical techniques such as finite volume and finite element methods; the trial and 
test spaces can be easily constructed, inhomogeneous boundary conditions and 
curved boundaries can be handled easily. They are also flexible in handling non¬ 
matching grids and in designing hp-adaptive mesh refinement. The dG methods 
were successfully applied for the steady state llT7ll^l48l . the time dependent lin¬ 
ear convection-diffusion-reaction |II1|2]], and the semi-linear steady state ll431l^ 
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optimal control problems. To solve the optimization problem, we here apply the 
optimize-discretize approach. The first order optimality conditions are derived 
and then they are discretized by using the dG method. We employ the projected 
nonlinear conjugate gradient (CG) method as an optimization algorithm limi . We 
show the controllability of the traveling waves of the FHN equation with target 
functions in the whole space-time domain and at the final time with and with¬ 
out sparse controls. In addition, the performance of the method is demonstrated 
for convection dominated problems by increasing the wave velocity, for sparse 
controls with different sparse parameters. 

Due to the semi-linearity of the FHN equation, the control problem is non- 
convex. Therefore, the fulfillment of the first necessary conditions does not im¬ 
ply the optimality. In order guarantee the optimality, the second-order sufficient 
optimality conditions (SSCs) have to be checked. In the recent years, the fulfill¬ 
ment of the SSCs for infinite dimensional and finite dimensional semi-linear PDF 
constrained optimal control problems has been investigated extensively (see, e.g., 
ifTlll for a recent overview). Except few examples with analytical solutions, it is 
not possible to prove the SSCs for the infinite dimensional problems since the un¬ 
known optimal solution is required. Therefore the finite dimensional approxima¬ 
tions of the infinite dimensional problem are considered. Provided that the local 
minima satisfies the SSC, one can check the SSC numerically by finding a bound 
for the distance between the local minima and discrete solution |l36l|33]. For this 
purpose, the associated coercivity constant of the reduced Hessian operator is esti¬ 
mated numerically by computing its smallest eigenvalue. Similar techniques were 
applied to measure how far the control obtained by a reduced order optimization 
model is away from a local full control solution ll^l26ll . Moreover, the Tikhonov 
regularization parameter in the cost function expresses the cost of the control, 
and it increases the numerical stability of the optimal solution. Recently the sec¬ 
ond optimality conditions have been investigated for semi-linear parabolic control 
problems with the objective function, not including the Tikhanov regularization 
term m. We test the discrete optimization problem for vanishing Tikhonov reg¬ 
ularization parameter as in ifTTll . The numerical results of the control of two di¬ 
mensional waves confirm the convergence of the optimal solutions for vanishing 
Tikhonov regularization parameter as it was demonstrated for the one dimensional 
wave solutions of the classical FHN equation in IfTTll^ . 

The paper is organized as follows: In the next section the optimal control 
problem governed by the convective FHN equation is described as a model prob¬ 
lem. We first prove the existence and uniqueness of a convective FHN equation, 
called as a state equation, by transforming into the one with monotone nonlin- 
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earity. Then, we introduce the first and second order optimality conditions. In 
Section we give a full discretization of the optimality system using the SIPG 
method in space and the backward Euler discretization in time. In Section]^ we 
discuss some benchmark examples with and without sparse controls. We investi¬ 
gate the effect of the Tikhonov parameter as it goes to zero. Further, with the help 
of the perturbation method and the smallest eigenvalue of the reduced Hessian, 
we find a bound for the distance between the local minima and discrete solution. 
The paper ends with some conclusions. 


2. Optimal control of the convective FHN system 

In this paper, we consider optimal control problems governed by the following 
convective FHN system: 


yt{x,t) - dyAy{x,t) + Y ■Vy{x,t) + g{y{x,t)) + z{x,t) = u{x,t) in Qt, (2.1a) 
Zt{x,t)-d^Az{x,t)+Y■'Vz{x,t) + e{z{x,t)-C 3 y{x,t)) = 0 in gr, (2.1b) 


with homogeneous Neumann boundary conditions 

d„y{x,t) = 0 , dnz{x,t) = 0 on Ej, 

Dirichlet boundary conditions 

y{x,t) =yD{x,t), z{x,t) = ZD{xd) onL^, 

and initial conditions 

y(v, 0 )=yoW, z{x,0)=zo{x) in Q. 


( 2 . 1 c) 

( 2 . 1 d) 

( 2 . 1 e) 


In this setting, let T > 0 be a fixed end time. We denote Qt the time space cylinder 
Qt :=£lx (0, T), where Q. = (0,L) x (0,if) is a bounded, Lipschitz domain. The 
lateral surface is denoted by E = T x (0,r). We use the notation E^ := Tz) x 
( 0 ,r) andE^ :=rjv X (0, T), where Dirichlet Tz) and Neumann boundaries, 
where the Dirichlet yD,ZD ^ and the Neumann boundary conditions are 

prescribed. Moreover, the initial functions are given as yo^zo G We denote 

the outward unit normal vector and the associated outward normal derivative on 
by n and dn, respectively. The diffusion coefficients are denoted by dy and d^. 
The parameters C 3 and e are real constants. Further, the function g(y) denotes the 
cubic polynomial nonlinearity 

^(y) =ciy(y-C 2 )(y-1) (2.2) 
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with the nonnegative real numbers c,-, / = 1,2, which is monostable for 0 < ci <20 
and C 2 = 0.02 ifT^ in contrast to the bistable cubic nonlinearity for the Schlogl 
equation |I7]|, the classical FHN equation ifTOll . and the diffusive FHN equation 
ll25]l . The velocity field denoted by V = {Vxi,Vx 2 ) is given along the xi-direction 
with a parabolic profile 

Vxi (X2) = aX2(H-X2), Vmax = a > 0, = 0, (2.3) 

where Vmax denotes the maximum wave speed of the velocity field. Moreover, the 
velocity field is the divergence free, i.e., div V = 0. 

We also make the following assumption for the solutions y, z 

0=:yo<y<yi, 0=:zo<z<zi a.e. 2r, (2.4) 

which is admissible for the sake of physical realism. We note that the bounds are 
constant. 

Here, we want to minimize an objective function of misfit type, i.e., the func¬ 
tion is designed to penalize deviations of the function values from the observed or 
measured data. We formulate our minimization functional such as 




(P^) min J{u) := I{u)+II j{u), 


(2.5) 


with 


I{u) = ^ j j {y{x,t)-yQ{x,tyf dxdt + ^ J J {z{x,t)-ZQ{x,t)f dxdt 


0 a 


0 Q. 


y z 

/ {y^^X)-yT{.xj)f dx+^ j {z{x,T)-zt{x,T))^ dx 

Q. n 

T 

j J {u{x,t))^ dx dt, 

0 Q. 


1 

;■(“) = II |M(x,t)| dx dt, 


0 Q. 


where the pair (y,z) denotes the solution of ( |2.1[ ) associated to control u. In 
(|2.1|), the partial differential equation for y is said activator equation, while the 
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one for z is called the inhibitor equation. The functions jtiZt G L°°(£2) and 
yQ-i^Q ^ are the given targets or desired states. We have given constants 

COq,(Oq,(Oj,(Oj, Tikfionov regularization parameter > 0, and sparse parameter 

jU > 0. 


We consider the optimization problem ( |2.5[ ) with pointwise box constraints 
u G 'Wad ■— {u £ L'^{Qt) ■ Ua < u{x,t) < U), for a.c {x,t) G Qt} (2.6) 


with the real numbers Ua<ui,. 

The aim of the optimal control is to ensure that the state variables j and z are 
as close as possible to the desired or observed states by minimizing the objective 
functional in the or Lfonorms. For well-defined optimal solutions, one has to 
show that there exists a unique solution (_y,z) of (2.1) for each u G Wad- The FHN 
equation ( |2.1| ) belongs to the class of semi-linear parabolic equations with a non¬ 
monotone nonlinearity. The theory of the existence and uniqueness of solutions 
(y,z) of the state equation ( |2.1[ ) is more delicate than the monotone nonlinearities 
lEiiia. Next section, we show the existence and uniqueness of the state equation 
( |2.1| ) by transforming the state equation to one with the monotone nonlinearity for 
the convective FHN equation p.l| ). 


2.1. Well-posedness of the state equation 

The existence and uniqueness of a weak solution for the Schlogel and FHN 
equation was shown in iQ by transforming ( |2.1| ) into the one with monotone non¬ 
linearity using the transformation y = e^^v with sufficiently large parameter r] 
0. Here we apply the same technique for the FHN equation with the convective 
term. Next, we construct the upper and lower solutions for the transformed equa¬ 
tion, which yield pointwise bounds for the desired solution following These 
bounds are then used as an initial iterates to construct two monotonically conver¬ 
gent sequences. Finally, we show that their common limit is the unique solution 
of the transformed equation with the monotone nonlinearity. 

Let us first perform the transformation of the state equation ( |2.1| ) by substitut¬ 
ing y = e^^v. Then, we obtain the following system: 

Vt — dyAv + Y •yv +e~^^g{e^^v) + rjv = e~^\u — z) in^r, (2.7a) 
Zt-dzAz-\-\-Yz + eiz-c-ie^^) = 0 mQj (2.7b) 

with the boundary conditions 

(9„v = 0, ^„z = 0 onEj, and e^*v = yD-, z = zd onE^, (2.7c) 
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and the initial conditions 


v(x,0) =};o(x), z(x,0)=zoW in (2.7d) 

Here, the nonlinear term 

g : V I—> e~^^g(e^^v) + rjv 

is a monotone non-decreasing function with respect to v for all (x, t) E Qj. More¬ 
over, it satisfies the following properties Sec. 4.3]: 

(i) For every fixed v e M is Lebesgue measurable in Qt. 

(ii) For almost all (x,t) E Qt, g twice continuously differentiable with respect 

to V and locally Lipschitz continuous of order 2 with respect to v, i.e., there 
exists L(p) = 6cie^^^ > 0 such that 

\gvv{x,t,Vi) - gvv{x,t,V2) \ <L{p)\vi -V2I 

holds with for all vi, V2 G M with |v, | < p, / = 1,2. 

The nonlinearity is uniformly bounded and monotone increasing in the following 
sense: 

(iii) There exists a constant C = ciC 2 + ri + 2e^^(ciC2 + ci) >0 such that 

|g(x, t , 0) I + 0) I + |gvv(^,^0) I < C. 


(iv) It holds 0 < gy(x,t,v) for almost all (x,t) E Qt, all v G M. 

Before defining a weak formulation of the system ( |2.7[ ), we need to define the 
following Hilbert space 

W{0J) ■.= {wEL^{<d,T\V)\ w EL^{0J■,¥*)} 


equipped with the norm 

IMk(o,T) = Q'dMOIlH IK Wi.)*)'. 


where V = and V* is the dual space of V. Now, we can define a weak 

solution of the system (|2.7|). 
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Definition 2.1. A pair of functions (v,z) G {W{0,T) nL°°(2r))^ is called weak 
solution of the system ( |2.7| ), if the equations 

j {vt,(p)v*,v dt + jj [dyVv-Vtp + y ■'Vv(p + gq)-e^\u-z)(p) dxdt =0, 

° (2.8a) 

J dt + jj V(p + V-Vj(p + £(2 — £-3f‘'''v)(p) dvd? =0, 

“ (2.8b) 


and 

v(-,0)=};o, z(-,0)=zo (2.8c) 

are satisfied for all (p G L^(0,r;if^(Q)). It is noted that V denotes the gradient 
with respect to x. 

Next, we give the definition of the ordered upper and lower solutions as done 
in (511321. The pair functions (v,z) and (v,z) in (1T(0, T) nL°°(27’))^ are said to 
be ordered upper and lower solutions of ( |2.71 , respectively, if they satisfy 

(v,z) < (v,z) 

and 

Vt — dyAv-\-\-Vv-\-g{x,t,v)—e~^\u — 'z)< 0 

z) — r/ 2 Az +V-Vz'+e(z —< 0 

dnV < 0 

dnZ< 0 

v< e-^^yo 

Z< ZD 

v{x,0)< yo(-^) 

z(x,0)< zo(-^) 

By taking 

v{x,t) = z(x, t) = M, v{x,t) = z{x, t) =0 


< Vt — dyAv + V • Vv 

<Zt- r/^Az + V • Vz 
+e(z-C3e^^v), 

< d„v, 

< dnZ, 

<V, 

<v{x,0), 

< z(v,0). 
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for some M > 0, we can rewrite our system ( |2.7| ) as 


Vf — JyAv + V -Vv + —e~^\u — M) 

Zt — + V- Vz + e(z — c^e^’M) 

dnV = 0 , dnZ 

V = e~^‘yD, z 
v{x,0) =yo{x), z(x,0) 


e-^\M-z), (2.9a) 

ecse^^M-v), (2.9b) 

0, (2.9c) 

ZD, (2.9d) 

zo(x)- (2.9e) 


Here, we have 


d{e-^‘{M-z)) 

dv 

d{e-^\M-z)) 

dz 


d{ec'ie^\M — v)) 
dz 

d(ec 3 e^‘ {M — v)) 
dv 


for all v,zE [0,M]. 

Now, we can state the existence and uniqueness of the system ( |2.7| ) for each 
control variable u. 


Theorem 2.2. Assume that the initial conditions yo <^nd zo are nonnegative func¬ 
tions, and ( |2.4| ) holds. Then, the system p.7| ) admits a unique solution (v,z) e 
(lT(0,r) nC^Qr))^ for each control u G ^ad- 

Proof We adopt the iteration technique introduced in ^37^ and construct sequences 
{(v*,f*)}^,^Q, {(y*,z with initial elements 

vO = v = M, z^=f=M, 
y_° = v = 0, z® = z = 0. 


Initiating from (y^^z^) and (y*,z^), [v^^^ and (y ^■'■^,z *+^) are found by 

solving 


V - dyAv + V • Vv + g{e^'v (m - M) 

z - d^Az ^+1 + V • Vz + e(z 

a„v*+i = o, a„z^+i 

v-^+i(x,0)=yo(x), z^+'(^,0) 


e-^fM-z'^), 

ec3e^\M-y^), 

0 , 

ZD, 

zo{x) 
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and 

vf+i-rf;yAv^+i+V-Vv^+i+g(en*+i)-e-^'(M-M) = e-^\M-z^), 
zf+i-J,Az^+i+V-Vz*+i+e(z*+i-C3e’^'M) = ec^e^^M-v^), 

a„v^+i=o, a„z*+i = 0, 

Y^+i=e"’^'}’D, Z *+1 = ZZ), 

v*+^(x,0) =};oW, z*+^(x,0) = zoW, 

respectively. The constructed sequence {(v^,z*)}^^q is monotone non-increasing 
and upper solution for all k. Conversely {(y z *) is monotone non-decreasing 

and lower solution for all k. Further, we have 

u^{x,t) < u^{x,t) and y*(x,r) < v^(x,r) 


for all A: G N and (x, t) £ Qj. 

By induction, we can verify the monotonicity of the sequence {v^}^=o' 
k = 0. 


VI - dyAv 1 + V • Vv 1 + g(e^^v ^)-e-^\u-M) 
z] — d^Az ^ -I- V• Vf ^ e(z ^ — c^e^^M) 

dnV ^ = 0, dnZ ^ 

v-i=e-^%, z-i 

v^(x,0)=yoW, 


ec^e^\M-y^), 

0 , 

zd, 

zoix). 


The property that v is an upper solution gives us 
V ° - dyAv 0 + V • Vv 0 + g{e^^v °) - {u - M) 

z0-rf^Az° + V-Vf0 + £(f0-c3e’^'M) 

5„V°>0, dnZ^ 

v^>e-^%, z-o 

v°(x,0)>yoW, z°(-^,0) 


> ec^e^\M-y^), 

> 0 , 

> Zd, 

> ZoW- 


Hence, we obtain 

v°-v/-4A(vO-vi)+V-V(vO-vi)+|(e’7fvO)-|(e’7^vi) > 0, 

z0-z,^-JzA(z^-z°)+V-V(z0-z1) + £(z0-z1) > 0, 

4(v-0_v-l)>0, dnif-Z^) > 0, 

(v^-v^)>0, (z°-zi) > 0, 

(v°-v ^)(x,0) > 0, (z°-z^)(x,0) > 0. 
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So it follows from the comparison principle for nonlinear parabolic equations v ° — 

V \ — z^ >0. Now if we assume that v*, z^~^ —z^ >0, one can easily 

show that V^ — V , z^ — z> 0. Analogously, the monocity of (y zcan 
be proved. 

Now, we show the convergence of the sequence (v z to a solution of (2.7). 
The sequence {v is monotone non-increasing and bounded from below by m = 
0. Hence by Lebesgue dominated convergence theorem |l5l |32]|, it converges to 

V in space LP{Q), p < oo. On the other hand, the sequence {z is monotone 
non-decreasing and bounded from above by z = M. It converges to z in a similar 
way. 

Finally, we prove the uniqueness of the solution of (2.7). Suppose that (vi,zi), 
(v 2 ,Z 2 ) are solutions of ( |2.8| ) and setv:=vi—V 2 , z:=zi—Z 2 - Then v and z satisfy 
the initial conditions obviously. Moreover, the following equations 


^ {vt,(p)v*,v dt+ {dyVv'V(p + Y-Vv(p + g{x,t,vi)(p 

— g(x,t,V 2 )(P — e^^(u — z)(p) dxdr = 0, 

[ {zt,(p)v*vdt+[[ (rfjVz-Vcp-l-V-Vzcp 
Jo JJqt 

+ e{z — C2e^^v)(p) dxdr = 0 


(2.10a) 


(2.10b) 


hold for all (p G VF(0, T)). Then, following the steps in ifTTll by taking (jp = v in 
and ^ = z in ( |2.10b[ ) we obtain the desired result v = 0 and z = 0. □ 

Hence, we can give the existence of an optimal control u for the optimal con¬ 
trol problem ( |2.5| ). 

Theorem 2.3. The optimal control problem ( |2.5| ) has at least one optimal solution 
u with associated optimal state y. 

Proof. Here we just sketch the key ideas of the proofs in iR^fTTl Sec. 5.3,Thm. 7.4] 
Since ^ad is non-empty and bounded in L°°{Qt), it is bounded in any space 
L^iQr) and it follows from the existence and uniqueness of the state variables 
that they are also bounded. Hence, the cost functional is bounded below, which 
allows the existence of an infimum. Therefore one can find a weakly convergent 
minimizing sequence due to the boundedness of this sequence. Then, the com¬ 
pact embedding results give us the strong convergence of the state in a weaker 
norm. Hence, there exists a feasible limit point, and convergence of the objective 
function can be shown using the continuity argument. □ 


(2.10a 
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Remark 2.4. In this paper, we show the well-posedness of the state equation by 
introducing upper and lower solutions for the state with monotone nonlinearity, 
obtained after a suitable transformation. However, there exists other possible 
reformulations in the literature, for instance, Schauder fixed point theorem applied 
in /IZ^/ for FitzHugh-Nagumo equation, Brouwer fixed point theorem applied in 
/f79l/ for a class of quasi-linear elliptic problems which are of nonmonotone type, 
Leray-Schauder fixed-point theorem applied in HWl for a coupled system of semi- 
linear parabolic reaction-diffusion equations, Faedo-Galerkin method applied in 
for a hyperbolic quasi-linear hemivariational inequalities. 

We continue this section by introducing the necessary and sufficient optimality 
conditions of the optimal control problem ( |2.5| ). 

2.2. First order optimality conditions 


The OCP 


min J{u) := f{yu,Zu,u) =I{u) + ilj{u) 


( 2 . 11 ) 




subject to the convective FHN equation ( |2.1| ) is a non-convex programming prob¬ 
lem so that different local minima might occur. Since any global solution is one 
of these local solutions, we set up the first-order optimality conditions satisfied by 
the local minima. 

The cost functional J(u) in ( |2. 11 ) consists of two terms with different smooth¬ 
ness. While the first part I{u) is smooth, the second part j{u) : L^{Q) —)• E is not 
differentiable, but it is subdifferentiable and the directional derivative is given by 



( 2 . 12 ) 


with 


dj{u) 



where 


?i{x,t) e 
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We note that the relations above are required only almost everywhere. With the 
help of the Lagrangian functional, we obtain the following variational inequality: 


/ (m)(v-m)> 0 


(2.13) 


that is, 

T 


x,t) + (Ouu{x,t) + flX(x,t)) (v(x,t) — u{x,t)) dxdt >0 Vv G 


0 Q. 


where p{x,t) with q{x,t) are called the adjoint variables as the solution of the 
following system 

-pt-dy^p-\■Vp + gy{y)p-ecj,q = ®Q(y-ye) ingr, (2.14a) 
-qt-dzl!^q-y -Vq + eq + p = (Oq{z-z^ in 2r, (2.14b) 

with the mixed boundary conditions 

dydnp{x,t) + {y■n)p{x,t) = 0, dzdnq{x,t) + (y■n)q{x,t) =0 onZ^, (2.14c) 
the Dirichlet boundary conditions 

p{x,t)=0, q{x,t)=0 onL^, (2.14d) 

and final time conditions 

p{x,T) = <xi^{y{x,T)-yTix)), q{x,T) = (0^[z{x,T) -zt{x)) inQ. (2.14e) 

The convection term in the adjoint system ( |2.14[ ) is the negative of the one in 
the FHN equation ( |2.1| ). As a consequence, errors in the solution can potentially 
propagate in both directions. Therefore, the numerical treatment of the state and 
adjoint systems together is more delicate. 

For CD„ > 0 and [i > 0, from the variational inequality ( |2.13| ) the following 
projection formulas are obtained lfmi38]l 

u{x,t) = |-^(/?(x,t)+juA(x,t))| fora.a (x,t) G 2r,(2.15) 

X{x,t) = P[-i,i] for a.a (x,t) G 2r, (2.16) 
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where the projection operator P[a/,] : M ^ is defined by 

P[a,/ 7 ](/) = max{a,min{/,&}}. 
Further, the following relation holds for almost all (x, t) G Qt 


u{x,t) = 0 


\p{x,t)\<ll, ifUa<0, 
p{x,t)>-ll, if Ua = 0. 


(2.17) 


We refer to |l9l [lOl |42|] for a further discussion on the projection operator and a 
derivation of (12.161) and p.lVI). 


2.3. Second order optimality conditions 

Due to the nonlinearity of the state equation ( |2.1| ), the optimization problem is 
non-convex. Therefore, the fulfillment of the first necessary conditions does not 
imply optimality. In order guarantee the optimality, the second-order sufficient 
optimality conditions (SSCs) have to be satisfied. The SSCs are related to certain 
critical cones that must be chosen as small as possible. 

Now, we give the critical cone related to our optimization problem for ©« > 0 
derived in ifTTIl : 

Cm = |v G I?‘{Qt) : V satisfies the sign condition and I (m)v-|- [i j (u) = o| 
with the sign condition: 

/ v{x,t) > 0, if u{xd) = 

\ v(x,t) < 0, if u{x,t) = U},. 

The set Cu is a convex and closed cone in I?{Qt). Moreover, if a is a local minima 
for (P^), then the following inequalities hold ifTTl Theorem 3.3], 

I (m)v^ > 0 VvgCm\{ 0}, equivalently / (m)v^ > 5||v||^2(q^) VvgCm, 5 > 0. 

(2.18) 

Then, under the assumption I (m)v^ >0 Vv G Cm\{ 0}, there exits 5 > 0 and 
ro > 0 such that 

J{u) + -\\v-u\\]2^^Q^^<J{y) Vm G t/arfnfiro(M), (2.19) 

where Br^iu) is the ball centered at u with radius tq. This shows the 

existence of a local minima, see HH Theorem 3.4] for details. 
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The verification of the SSCs is difficult because the solution of the infinite 
dimensional problem is required. Even if it is known, it would still be tedious to 
check that the SSC holds, since it requires the exact solution of linearized PDEs. 
However, there exist some numerical studies on the SSCs, see e.g., Il^l29ll^ . 
Here we determine the constant 5 by computing the smallest eigenvalue of the 
reduced Hessian as introduced in ll29l of the finite dimensional dG discretized 
OCR 

Now, we can state the following theorem to measure the distance between 
the local minima u and the local discrete solution m/, obtained by applying the 
discontinuous Galerkin discretization for spatial discretization and the backward 
Euler for temporal discretization in Section]^ 


Theorem 2.5. Let ube a local minima 0/(0. Assume that u satisfies the second- 
order sufficient condition ( |2.i<Sp and ( |2.i9p . Ifuh is the discrete solution such that 
Uh G Brfiu), then it holds 




( 2 . 20 ) 


where a perturbation function is defined as the following 


CW := 


and 0 < 8 < 8. 


- min{0, OiuUh + p/, + pA/,}, 

~ + PhA- pA/i), 

- max{0, (OuUh + Ph + ftA/,}, 


If ~ 
f ^ 
f 


Proof Let m/, be a discrete solution that need not be optimal for the continuous 
problem (2.5), and let ph and Xh be the associated adjoint and sparse. If Uh were 
optimal, then ph + + juA/, = 0 should be satisfied in almost all points x G H, 

where Ua < u^ < uy holds. If not, then ph + (Ou^h + pA/, + ^ = 0, where ^ is a 
perturbation function, adopted from ||3]|. Although a/, is possibly not be optimal 
for ( |2.5| ), it is optimal for the perturbed optimization problem 


min 7(M) + (C,M)L 2 (g^ 

^ad 


Inserting u in the discrete variational form and u^ in the continuous discrete vari¬ 
ational form, we obtain 


(m) + — wj > 0. 


(2.21a) 

(2.21b) 
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Addition these equations in (|2.21|) gives us 


(j'(Uh) - j'{u),u - Uhj + {^,U-Uh) > 0. 


By the mean value theorem, we obtain 

— J {u){u —Uh) + {^,U —Uh) >0 


( 2 . 22 ) 


(2.23) 


for some m e {v G Uad ■ v = u-\-t{uh — «), t G (0,1)}. Then, when we apply the 


SSC (2.18) with Cauchy-Schwarz inequality, we find 


2 

5 ||m/i — < IICIIl2(2j-)II“/! ~ 

which is the desired result. 


(2.24) 

□ 


Here we follow ll^ Remark 3.3] and select 5' := 5/2. If «/, belongs to the 
neighborhood of u, we can estimate in the following 


'-Uh\y{QT) < ^IICIlL2(er)- 


(2.25) 


Remark 2.6. We also remark that the presence of the so-called Tikhonov param¬ 
eter ^ in the cost functional {2.5) is extremely important. The standard second- 


order optimality conditions do not hold for vanishing Tikhonov parameter (Ou = 0. 
By introducing new critical cones, the second-order optimality conditions are es¬ 
tablished in 4771/ for sparse optimal control governed by FitzHugh-Nagumo sys¬ 
tem, and in 4771/ for general nonlinear functions. Hence, the new second order 
conditions are used for proving the stability of locally optimal solutions with re¬ 
spect to (Ou 0. This theoretical result was confirmed for one dimensional wave 
solutions of the classical FHN equation with sparse controls in 47711771/ . In this 
paper, we apply the theory introduced in 4771 1771/ for two dimensional wave solu¬ 
tions of the convective FHN equations. 


3. Discretization in space and time 


We discretize the state system ( |2.1| ) and the adjoint system ( |2.14| ) by the sym¬ 
metric interior penalty Galerkin (SIPG) in space and the backward Euler dis¬ 
cretization in time. 
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3.1. SIPG discretization of the state and adjoint equations in space 

The interior penalty dG methods are well established in the literature and the 
details can be found in the classical texts like Bl l^ . 

We assume that the domain Q is polygonal domain. We denote {ffijh as a 
family of shape-regular simplicial triangulations of Q., see, e.g., ifTdll . Each mesh 
Ifh consists of closed triangles such that Q. = ^ holds. We assume that the 

mesh is regular in the following sense: for different triangles e 

the intersection Ki n Kj is either empty or a vertex or an edge, i.e., hanging nodes 
are not allowed. The diameter of an element K and the length of an edge E are 
denoted by Hk and hs, respectively. 

We split the set of all edges #/i into the set of interior edges, the set 
of Dirichlet boundary edges and the set of Neumann boundary edges so that 
= <^/f U with U . Let n denote the unit outward normal to dO.. 

For the activator y and the inhibitor z, we define the inflow and outflow parts of 
do. by 

r- = {x G : V(x) • n(x) < 0} , r+ = {x G : V(x) • n(x) > 0} . 

Similarly, the inflow and outflow boundaries of an element K are defined by 

dR- = {xedK\ V(x) ■ ni^(x) < 0}, = {xedK: V(x) • (x) > 0}, 

where n/f is the unit normal vector on the boundary dK of an element K. 

Let the edge E be a common edge for two elements K and K‘^. For a piecewise 
continuous scalar function y, there are two traces of y along E, denoted by y|£: 
from inside K and y^\E from inside . Then, the jump and average of y across 
the edge E are defined by: 

[tv]] {{3^}} = 

Similarly, for a piecewise continuous vector field Vy, the jump and average across 
an edge E are given by 


For a boundary edge E G ^£2, we set {{y}} = y and [[y]] = yn. 

Recall that for dG methods, we do not impose continuity constraints on the 
trial and test functions across the element interfaces. As a consequence, the 
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weak formulation must include jump terms across interfaces, and typically penalty 
terms are added to control the jump terms. Then, we define the discontinuous dis¬ 
crete space as follows: 

Wh = {we l 2(Q) : w l^e F\K), yK e =^4 , (3.3) 


where (K) is the set of piecewise linear polynomials defined on K. We note that 
the space of discrete states and the space of test functions are identical due to the 
weak treatment of boundary conditions for dG methods. 

Then, the semi-discrete formulation of the state system ( |2.1| ) for Vw e Wh and 
t G (0, T] becomes 


dyh 

dt 


,w]yah,y{yh^^) + h,y{yhM + Ch,z{zh,w) = 4,3,(w) + (m/,, w),(3.4a) 


= {yo,w), 

+ ah,z{zh,w)+ bh,z{Zh , w) -f Ch,y{yh , w) = £h,z(w), 

{zhi-, 0 ),w) = (zo,w)> 


(3.4b) 

(3.4c) 

(3.4d) 


where the (bi)-linear terms are defined for i = y,z and Vw G W), 

CltA'’,'"’) = Y, jdyv-Vwdx: 

(j,-vv}}H + {W-vm.}}.h)<;s 

ds+ ^ Iw-Vvwdx 

Ke^hi 

+ E / V-n(v^-v)wJ 5 - £ / \-nvwds, ( 3 . 5 a) 

4 ,/(>v) = Y, j • M “ ds 

- Y j y-^iDwds, ( 3 . 5 b) 


- I 

+ I ^/h-H 


Ke£^h 


dK-nr- 
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(3.5c) 



Ch 


Ch. 


’{yM= L -£c3ywdx, 

ziz,w)= Y, Jzwdx, (3.5d) 


where the parameter <y G Mq is called the penalty parameter which should be 
sufficiently large to ensure the stability of the dG discretization; independent of 
the mesh size h and of the diffusion coefficients d, as described in ll35l Sec. 2.7.1] 
with a lower bound depending only on the polynomial degree. Large penalty 
parameters decrease the jumps across element interfaces, which can affect the 
numerical approximation. However, the dG approximation can converge to the 
continuous Galerkin approximation as the penalty parameter goes to infinity (see, 
e-g-, imi for details). 

For each time step, we can expand the discrete solutions of the states y,z, and 
the control u as 


N ri]. 


N ri]. 


(0 = L L yj^j ’ zh{t) = Y'Y ’ 


/=l7=l 


/=l7=l 


N rik ^ 

and Uh{t) = YH 
i=lj=l 


where yj, Zj, Uj, and are the unknown coefficients and the basis functions, 
respectively, for y = 1,2, • • • ,% and i = 1,2, • • • ,A^. The number N denotes the 
number of dG elements and % = (k+ l)(k + 2)/2 is the local dimension of each 
dG element depending on the order k of the polynomial basis. 

Inserting p.6| ) into p.4| ), we obtain 

M—+ A^y + by (y) + = ly + Mw, (3.7a) 

M^+A^z + B^z + Cyy = 1„ (3.7b) 

where y, z and u are the unknown coefficient vectors, i.e.. 


z 

u 


= 


yuk 


5 5 


nj,) 
Z^) 


_ ^N\ 

5 ■ ■ ‘ 5 ’ ’ ’ ’ ’ ’ ’ ’ ’ ^rik)^ 


rik) 


M is the mass matrix. Ay and A^ are the stiffness matrices corresponding to 
(^h{yhiw) and ah{zhiw), respectively, and ly and are vectors corresponding to 
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£/i 3 ,(w) and £h^^{w), respectively. Moreover, Cy, and are the matrices cor¬ 
responding to Ch^y, bh^^, and Ch^^, respectively. Further, by(y) is the vector corre¬ 
sponding to the nonlinear part bh,y{y,w). 


In the literature, there exist two approaches for the solution the OCP ( |2.5| ) nu¬ 
merically, i.e., discretize-then-optimize (DO) and optimize-then-discretize (OD). 
In the DO approach, the objective function ( |2.5[ ) and the state equation ( |2.1| ) are 
discretized first, and then the discrete optimality system is formed. In the OD 
approach, the optimality conditions consisting of the state system CT. the ad¬ 
joint system ( |2.14[ ) and the variational form ( |2.13[ ) are derived. Then, the infinite 
dimensional optimality system is discretized. We here apply the optimize-then- 
discretize approach. The discretization and optimization commute for the SIPG 
discretized linear steady state convection-diffusion-reaction OCPs |l47l|48]|. For 
linear time dependent OPCs the OD and DO approaches commute for the dG dis¬ 
cretization in time and the SIPG discretization in space Q. The commutativity 
property is no more valid for semi-linear elliptic and parabolic OCPs. 

In the optimize-then-discretize approach, we derive directly semi-discrete form 
of the adjoint equations ( |2.14| ) using the SIPG discretization as for the state equa¬ 
tions ( |3.4| ): 

dph 


dt 

dqh 


dt 


+ah,p{ph,w)+bh,p{ph,w)+Ch,q{qh,w) +ih,p{w), 

{Ph{-,T),w) =(ol {{yh{-J)-yT{x)),w ), 

+ah,q{.qh,w)+bh,q{qh,w)+Ch,p{ph,w) +£h,q(w), 

{qhi-,T),w) {{zh{-,T)-ZTix)),w ), 


where the bilinear forms a/, p and a^ q are similar to the state ones (|3.5|) with the 
negative convection terms, i.e., — V(x,y). As an extra term, they just contain the 
contribution of the mixed boundary conditions, i.e.. 


/ {Y-n)phwds and / {Y-n)qhwds. 

Ee<^hE 


The other terms are given by 


bh,piP^^)= E 8y{y)pwdx, Ch,q{q,w)= -£C 3 qwdx, (3.9a) 
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bh,qiq,w)= Y, Ch,p{p,w)= Y pwdx, (3.9b) 

Ke-^hK K(^^hi 

4,p(>v) =-0^ jQwdx, =-Q)^ Y ZQwdx. (3.9c) 

Ke^hi Ke^hi 

3.2. Time Discretization 

We solve the nonlinear system of ordinary differential equations ( |3.7| ) in time 
by the backward Euler method. We first divide the time interval [0,7] 0 = to < 
t\ < ... < t^j = T, with step size tn = tn — tn-i, n = 1,2 ,...,Nt and then obtain 
the following fully discretized state system: 

—— y„_i) + + by(y„) + + (3.10a) 

'^n 

]Vl(Zn Zfi—l) ”1” AjZn + (3.10b) 

'^n 

for n= 1,2,... ,Nt. We solve the nonlinear system by Newton’s method. 

We discretize the adjoint variables p and q analogously to the state variables 
y, and z as 


Nnj,^ Nrik 

Ph{t) = YYPj^j ^h{t) = YY^j^j- 

/=l7=l i=l7=l 

By inserting p.ll| ) into ( |3.8[ ) and applying the backward Euler discretization in 
time, we obtain the following fully discretized adjoint system: 

Pn^ -1~^pPn—i 3“Bp(yn)pn —1 =n)Q(My^-f , (3.12a) 

'^n 

M(^n—1 ^qqn—1 3~^qqn—l 3~ CpPfi—i =tOQ (MZji 3“ 1^) (3.12b) 

'^n 


for n = Nt, ... ,2,1. 
forms ah,p and ah,q 


The matrices A.p,Aq correspond to the discretized bilinear 
respectively, whereas the matrices Bp(y), Cp, B^, Cq and the 


vectors Ip,l^ are the discretized forms in (3.9), respectively. 


4. Numerical results 


In this section, we provide numerical results for the OCP governed by the 
convective EHN equation ( |2.5| ). There exists several optimization algorithms for 
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OCPs governed by semi-linear equations. We have used in this paper the projected 
nonlinear conjugate gradient (CG) method EOl l2T]l . which was applied to the 
Schlogl and FHN equations ||71[l0l[II]|. The projected nonlinear CG algorithm is 
outlined below: 


Initialization Choose an initial guess uq, an initial step size and stopping tol¬ 
erances Toll and T0I2. Then, compute 




initial states (yo,zo) = t>y solving ( |3.10| ), 

initial adjoints {po,qo) = (Pyo,zo>^>'o,zo) solving ( |3.12| ), 
subgradient of j by ( |2.16[ ), i.e., Xq = Kq.po, 
subgradient of J, i.e, go = + Po + 

anti-subgradient of J, i.e., do = —go. 


Set k:=0. 


Step 1. (New subgradients) Update 




control, i.e., = Uk + Skdk, 


states {yk+i,Zk+i) = {yuk+i,Zuk+i) by solving ( |3.10| ), 
adjoints (pk+uqk+i) = {Py,+uZk+i^^yk+uZk+i) solving ( |3.12| ), 
subgradient of j by (2.16), i.e., Xk+\ = A, 


'Uk+\^Pk+\ ’ 


subgradient of J, i.e, g^t+i = (OcUk+i + Pk+\ + ^h+\- 


Step 2. (Stopping Criteria) Stop if ||gy(;+i|| < Toll or \\Jk+\ —JkW < T0I2. 

Step 3. (Direction of descent) Compute the conjugate direction fik+i according 
to one of the update formulas such as Hestenes-Stiefel, Polak-Ribiere, Fletcher- 
Reeves, and Hager-Zhang, see e.g., |l7l|2T]] for details. 


dk-\-\ = -gk+\ -\-^k+\dk- 


Step 4. (Step size) Select step size ^jt+i according to some standard options such 
as bisection, strong Wolfed-Powell, see e.g., Q |20]| for details. Set k =: 
k-\- \ and go to Step 1. 

For the computation of the reduced Hessian, we use the BFGS algorithm ll3Tl 

El: 
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• Set Hq = 1. 

• Update for A: = 1,2,... 


Hk+\ 


^kqj 

qlrk 


{Hkrk){HkrkY 

rlHkSk 


where = Wfe+i - uj, and = gk+\-gk- 


After setting the reduced Hessian matrix at each time level, we compute the small¬ 
est eigenvalue of the Hessian by using MATLAB® eig. 

In all numerical experiments, we started with the initial control m = 0. The 
optimization algorithm is terminated, when the gradient is less than Toll = 10“^ 
or the difference between the successive cost functionals is less than T0I2 = 10“^. 
We use uniform step sizes in time At = 0.05. The penalty parameter o in p.5[ ) is 
chosen as (7 = 6 on the interior edges and cr = 12 on the boundary edges. 

If it is not specified, we use the following parameters in all numerical examples 


Cl = 9, C 2 = 0.02, C 3 = 5, e = 0.1, Ji=c ?2 = l 

on a rectangular box = [0, L] x [0, if] with L = 100 and ii = 5 as in ifTbll . Further, 
the final time is taken as T = 1. 

All simulations are generated on Intel(R) Core(TM) i7-4720HQ CPU@2.60GHz, 
16GB RAM, Windows 8, with MATLAB® R2014a (64-bit). 


4.1. Optimal control in the space-time domain 

We first consider the OCP with the desired state functions defined in the whole 
space-time domain Qj with ©g = cOg = 1, (Oj = (Oj = Q, and ©„ = 10“^. The 
step sizes in space are taken as Axi = Ax2 = 0.5. The desired states are chosen as 
the solution of uncontrolled FHN equation (|2.1|) 


^ [0, otherwise, 'e.y 1 y otherwise. 


with the initial conditions given by 


f 0.1, if 0<xi <0.1, 
[ 0, otherwise. 


zo(-^,0) =0. 


Further, we define the admissible set of controls as 


^ad •= {u G L°°{Q)' —0.2 < u{x,t) < 0 for a.e {x,t) G Q}. 
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Figure 1: Example [43| optimal controls u (top) and associated states (bottom) at r = 0.75 without 
sparse control for ymax=16, 32, 64, 128 (from left to right). 


^max 

J 

#ite. 

#search 

#Newton 

16 

2.91e-2 

64 

245 

787 

32 

5.80e-2 

77 

297 

956 

64 

1.16e-l 

96 

373 

1203 

128 

2.30e-l 

123 

481 

1554 


Table 1: Example [43] cost functional /, number of nonlinear CG iterations, line searches, Newton 
steps without sparse control. 




Figure 2: Example 4.1 profiles of optimal controls u and associated states along xi axis at 


t = 0.15 without sparse control for various values of Vm 
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(a) Control u 


(b) State y 


4.1 


profiles of optimal controls u and associated states >> along xi axis for 


Figure 3: Example 
Vmax = 64 without sparse control for various values of ci 



Figure 4: Example |4.l[ optimal controls u (top) and associated states y (bottom) for the sparse 
parameters ^ = 1/500,1/100,1/50, 1 /35 (from left to right) and Emax = 32 at f = 0.75. 

We first investigate the numerical solutions of the optimization problem ( |2.5| ) 
without the sparse parameter, i.e., /t = 0. Figure [T] demonstrates the computed 
solutions of the control u and their associated states _y at t = 0.75 for various 
values of the Vmax, respectively. The control u, bounded by the box constraints, 
exhibits the same behavior of the state y. When the value of Fmax increases, the 
controlled solutions become more curved, as for the uncontrolled solutions ifT^ . 
Optimal values of cost functional J, and the number of iterations, line searches, 
and Newton steps in Table [T] are increasing for the higher values of Vmax because 
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(a) Control u (b) State y 


Figure 5: Example 4.1 


profiles of optimal controls u and associated states y at t = 0.75 with 


Vmax = 32 for various values of /i. 


the linear system of equations to be solved become more stiff due to convection 
dominated character of the OCP problem. 



J 

At; 

#ite. 

#search 

#Newton 

1/35 

l.lOe-0 

1.30e-l 

170 

666 

2193 

1/50 

7.31e-l 

8.21e-2 

156 

610 

1971 

1/100 

3.60e-l 

3.47e-2 

121 

470 

1516 

1/500 

1.21e-l 

5.93e-3 

88 

341 

1099 

0 

5.80e-2 

0 

77 

297 

956 


Table 2: Example [4l| optimal values of J, jx j, and number of iterations, line searches, and Newton 
steps for Vmwc = 32. 


One of the important features of the optimal control is the robustness of the 
control with respect to the parameters. However, in the PDE-constrained opti¬ 
mization context, we do not always have an explicit control function as in ll34ll^ . 
In that case, we check the robustness of the control by solving the optimal control 
problem for different values of the parameters. For this example, we study the 
effect of the maximum velocity Vmax on the controlled wave solutions and of the 
parameter ci on the stability of the waves. As V^ax increases, the wave profiles are 
elongated xi directions and control is adjusted accordingly by shifting to the right 
in Figure]^ Figure|^shows that the solutions display oscillations for higher values 
of Cl as the uncontrolled solutions in ifT^ . Since the control variable stays inside 
the bounds of the control, we can conclude that the control variable is robust with 
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respect to variations of the system parameters. 

Now, we look into the effect of the sparse parameter ji on the optimization 
problem ( |2.5| ). Higher values of jl cause the sparsity of the optimal control, which 
can be clearly seen from Figures]^ and Figure]^ shows profiles of the controls 
u and associated states y along x\ direction for various values of jU. The values 
of the optimal cost functional and number of iterations with line searches and 
Newton steps are given in Table When the sparse parameter ji increases, the 
smooth and non-smooth cost functionals increase as well as the number nonlinear 
CG iterations, line searches and Newton iterations. However, we will show in the 
next example with terminal control, the number of nonlinear CG iterations, line 
searches and Newton iterations can decrease when the sparse parameter increases. 
Similar behavior was observed in ifTOll for different problems with sparse controls. 




(a)p = 0 (b)p = 1/100 


Figure 6: Example 4.1 L^{Q) errors for -Jrefll’ W^oh “^refll, and -Mref||. 


Finally, in this example we study the convergence of the optimal control and 
its associated state for vanishing Tikhonov parameter ©„ for ;U = 0 and ft 0. 
The convergence of optimal states and control for vanishing Tikhonov parameter 
(Ou plays an important role for checking the SSCs according to the theory in ifTTl 
[13]]. It was demonstrated for the one dimensional wave solutions of the classical 
FHN equation with sparse controls in liTTl 1^ . Here we obtain similar results 
for two dimensional wave solutions of the convective FHN equations as shown in 
Figure]^ 

We fix the maximum velocity Vmax = 64. We take y^gf := yie_io, Zref := zie-io. 
and Mref := Mie -10 as reference solutions to determine the order of convergence as 
(Ou i 0. We look for the errors, i.e., -yrefll’ ll^co^ -Zrefll. and - Wrefll as 
(Ou i 0. We observe that the errors in L^-norm decay linearly as (OuiO in Figure]^ 
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0 

fi = 1/100 

(Ou 

Wycou-yaWma) 


\Zo)u ~Zq\\l2(Q) 

bcou-yaWL^Q) 

lk®c ^q\\l'^{q) 

1 

5.77167e-l 

1.09495e-l 

5.28730e-l 

9.90557e-2 

le-1 

3.84021e-l 

7.32190e-2 

4.48344e-l 

8.48398e-2 

le-2 

3.42134e-l 

6.53493e-2 

4.02207e-l 

7.44464e-2 

le-3 

3.39051e-l 

6.46970e-2 

3.95422e-l 

7.25312e-2 

le-4 

3.38605e-l 

6.46413e-2 

3.94826e-l 

7.23321e-2 

le-5 

3.38590e-l 

6.46329e-2 

3.94671e-l 

7.22975e-2 

le — 6 

3.38588e-l 

6.46321e-2 

3.94668e-l 

7.22959e-2 

le-7 

3.38588e-l 

6.46320e-2 

3.94667e-l 

7.22957e-2 

le-8 

3.38588e-l 

6.46320e-2 

3.94667e-l 

7.22957e-2 

le-9 

3.38588e-l 

6.46320e-2 

3.94667e-l 

7.22957e-2 

le-10 

3.38588e-l 

6.46320e-2 

3.94667e-l 

7.22957e-2 


Table 3: Example 4.1 L?-{Q) errors for - Jel| and Wzo)^ -zq\\. 


for both cases jU = 0 and jJ, = 100. For the norm of L°°{Qt) this looks similar. 

The trajectories of the target functions, i.e., yQ and zq, are quite achieved as shown 
in Table E 

4.2. Terminal control 

Our second test example is that the target functions, that is, yji-^T) and 
are only given at the final time. Regularization parameters are chosen 
d& oyQ = (Oq = Q, (Oj = (Oj = \, and ©„ = 10“^. We take the step sizes in space 
as Axi = Ax 2 = 0.125. We construct the desired functions as done in the previous 
example. They are given by 


yT{x,T)=yr,aA{x,T/I) and zt{x,T) = z^vax^x,!/I), 

where ynat and Znat are the solutions of the uncontrolled convective FHN equation 
at the final time T = \ with the initial conditions 


where Xa = 2 and xt = 2.2. The admissible set of controls is defined as 
•= e L°°{Q) '■ 0 < u{x,t) < 0.2 for a.e {x,t) G Q}. 
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Figure 7: Example |43| optimal controls u (top) and associated states y (bottom) at f = 1 without 
sparse control for V'max=16, 32, 64, 128 (from left to right). 


^max 

J 

#ite. 

#search 

#Newton 

16 

1.12e-2 

15 

61 

198 

32 

4.30e-2 

16 

65 

211 

64 

1.42e-3 

22 

89 

289 

128 

2.89e-4 

25 

95 

304 


Table 4: Example [43| optimal value of cost functional J, and number of nonlinear CG iterations, 
line searches, and Newton steps without sparse control. 


We observe similar behavior for the controls and associated states in Figure]^ 
and for the optimal value of cost functional J, the number of nonlinear CG it¬ 
erations, line searches, and Newton steps in Table as obtained for the whole 
space-time domain Example ( |4.1| ). 

The effect of the sparse parameter for the controlled solutions is displayed 
in Figure and As stated before, in this case for higher values of the sparse 
parameter ji the number of nonlinear CG iterations, line searches and Newton 
steps decrease as given in Table We observe the same linear decay for the states 
and the control for vanishing Tikhonov parameter in Figure [T^ as in the whole 
space-time cylinder, cf. Figure]^ The values of — yr|| and \\zo)u ~ zt\\ in L^- 
norm are given in Table We observe that the features of the desired trajectories 
are quite achieved. 
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Figure 8: Example |4. 2 [ optimal controls u (top) and associated states >> (bottom) at f = 1 for the 
sparse parameters = 1/2000,1/200,1 /150, and 1/100 (from left to right) and Fmax = 64. 



(a) Control u 


Figure 9: Example [4.2| profiles of optimal controls 
various values of the sparse parameter fi. 



u and associated states y (right) at r = 1 for 


Finally, in this example, we check the distance of discrete control a/, from the 
local minima u for the verification of SSCs ( |2.19| ). Since we do not know the exact 
optimal solutions of this example, we take the discrete solutions computed with 
Ax\ = Ax 2 = 0.3125 as reference solutions. Table |7] shows the numerical errors 
of ||m — M/,11 and the error estimates given in ( 2.25| ) for ft = 1 /200. The numerical 
results verify the Theorem |2.5[ However, the error and estimator are not reduced 
sufficiently for finer discretizations. 
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Jopt 

f^jopt 

#ite. 

#search 

#Newton 

1/100 

1.37e-l 

1.15e-2 

2 

5 

16 

1/150 

9.24e-2 

7.70e-3 

3 

7 

18 

1/200 

7.00e-2 

5.88e-3 

4 

8 

19 

1/2000 

8.79e-3 

6.43e-4 

19 

74 

238 

0 

1.42e-3 

0 

22 

89 

289 


Table 5: Example 4.2 optimal value of /, jlj, and number of nonlinear CG iterations, line 


searches, and Newton steps for V^ax = 64. 




(a)At = 0 (b) At = 1/200 


Figure 10: Example 


4.2 


HG) errors for \\y^^ - Jrefll^ W^ok, -Zrefll, and ||Ma,„ - Wrefll- 


5. Conclusions 

Numerical results of the optimal control governed by the convective FHN 
model with traveling waves reveal different aspects of the parabolic semi-linear 
optimal control problems investigated. The second order optimality conditions 
for local solutions in form of 2D traveling waves are verified numerically for van¬ 
ishing Tikhonov regularization parameter as done for one dimensional waves 
of the classical FHN equation in ifTTlI^ . By using the second order optimality 
conditions, we estimate the measure between the discrete solution and the local 
minima. The control of the convection dominated problems with wave solutions 
requires a large amount of computing time. In a future study we will investi¬ 
gate the reduced order modeling with proper orthogonal decomposition (POD) in 
space and model predictive control (MPC) in time |l39]|. 


31 


























M = 


0 

jU = 1/200 

(Ou 

\\yco,-yT\\LHa) 



lb®„-3^r||L2(n) 

\\Z(Oc-Zt\\l'^(q,) 

1 

6.53791e-2 

2.77918e-2 

7.92435e-2 

2.81181e-2 

le-1 

5.21169e-2 

2.74520e-2 

7.68586e-2 

2.80154e-2 

le-2 

4.87670e-2 

2.73638e-2 

7.13606e-2 

2.7785 le-2 

le-3 

4.86308e-2 

2.73684e-2 

7.11792e-2 

2.77772e-2 

le-4 

4.86178e-2 

2.73689e-2 

7.11610e-2 

2.77764e-2 

le-5 

4.86165e-2 

2.73690e-2 

7.11591e-2 

2.77763e-2 

le — 6 

4.86164e-2 

2.73690e-2 

7.11589e-2 

2.77763e-2 

le-7 

4.86164e-2 

2.73690e-2 

7.11589e-2 

2.77763e-2 

le-8 

4.86164e-2 

2.73690e-2 

7.11589e-2 

2.77763e-2 

le-9 

4.86164e-2 

2.73690e-2 

7.11589e-2 

2.77763e-2 

le-10 

4.86164e-2 

2.73690e-2 

7.11589e-2 

2.77763e-2 


Table 6: Example 4.2 IJ-(Q) errors for ||jg, — j 7 -|| and \\z(Ou 


Axi = Ax2 

\\U-Uh\\ 

IIICII 

2.5 

2.43e-2 

5.92e-2 

1.25 

2.43e-2 

5.60e-2 

0.625 

2.46e-2 

4.78e-2 


Table 7: Example |4.2[ Numerical errors of 11 m — m/, 11 and error estimates § 1 1 Cl I for sparse controls 
with /t = 1 /200. 
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